/*
 * Copyright 2009-2019 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
#define BOOST_TEST_MAIN

#define BOOST_TEST_MODULE gw_test
#include <boost/test/unit_test.hpp>
#include <fstream>
#include <votca/xtp/gw.h>

using namespace votca::xtp;
using namespace std;

BOOST_AUTO_TEST_SUITE(gw_test)

BOOST_AUTO_TEST_CASE(gw_full) {

  ofstream xyzfile("molecule.xyz");
  xyzfile << " 5" << endl;
  xyzfile << " methane" << endl;
  xyzfile << " C            .000000     .000000     .000000" << endl;
  xyzfile << " H            .629118     .629118     .629118" << endl;
  xyzfile << " H           -.629118    -.629118     .629118" << endl;
  xyzfile << " H            .629118    -.629118    -.629118" << endl;
  xyzfile << " H           -.629118     .629118    -.629118" << endl;
  xyzfile.close();

  ofstream basisfile("3-21G.xml");
  basisfile << "<basis name=\"3-21G\">" << endl;
  basisfile << "  <element name=\"H\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"5.447178e+00\">" << endl;
  basisfile << "        <contractions factor=\"1.562850e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"8.245470e-01\">" << endl;
  basisfile << "        <contractions factor=\"9.046910e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"1.831920e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "  <element name=\"C\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"1.722560e+02\">" << endl;
  basisfile << "        <contractions factor=\"6.176690e-02\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"2.591090e+01\">" << endl;
  basisfile << "        <contractions factor=\"3.587940e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"5.533350e+00\">" << endl;
  basisfile << "        <contractions factor=\"7.007130e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << endl;
  basisfile << "      <constant decay=\"3.664980e+00\">" << endl;
  basisfile << "        <contractions factor=\"-3.958970e-01\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"2.364600e-01\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"7.705450e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.215840e+00\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"8.606190e-01\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << endl;
  basisfile << "      <constant decay=\"1.958570e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "</basis>" << endl;
  basisfile.close();

  Eigen::VectorXd mo_eigenvalues = Eigen::VectorXd::Zero(17);
  mo_eigenvalues << -10.6784, -0.746424, -0.394948, -0.394948, -0.394948,
      0.165212, 0.227713, 0.227713, 0.227713, 0.763971, 0.763971, 0.763971,
      1.05054, 1.13372, 1.13372, 1.13372, 1.72964;
  Eigen::MatrixXd mo_eigenvectors = Eigen::MatrixXd::Zero(17, 17);
  mo_eigenvectors << -0.990744, 0.205494, 1.12419e-14, 1.45929e-14, 3.75871e-15,
      0.169799, 8.81326e-15, -7.59696e-14, 4.96934e-12, 3.56824e-13,
      -4.00071e-15, -9.34192e-15, -0.116538, -1.31882e-14, 3.3298e-14,
      -4.67337e-12, 0.0231349, -0.089086, -0.203519, 1.36072e-14, 4.0936e-14,
      -2.48734e-12, -0.0756531, 8.87485e-15, 1.17892e-14, -6.66069e-12,
      -3.61695e-12, -5.11501e-14, -5.65398e-14, 0.160728, 4.94743e-15,
      8.66668e-14, -1.06301e-11, 1.99744, -1.87806e-15, -2.96718e-13,
      -0.0253717, -0.29342, 0.203621, 1.37293e-12, 0.0616491, 0.221397,
      -0.161522, 0.405432, 0.0624782, 0.566771, -2.44075e-11, -0.205261,
      -0.709949, 0.521155, 5.26171e-12, -1.61071e-14, -4.43669e-13, 0.263897,
      0.122148, 0.208899, 1.16501e-12, 0.160859, -0.163271, -0.162399, 0.402943,
      -0.523447, -0.230538, -2.4163e-11, -0.512016, 0.531592, 0.522505,
      5.1261e-12, -8.74496e-15, -3.64573e-13, -0.240655, 0.164879, 0.207606,
      1.26241e-12, -0.221879, -0.0568539, -0.162616, 0.403447, 0.460008,
      -0.33931, -2.43511e-11, 0.716567, 0.176478, 0.522636, 5.19033e-12,
      0.0813118, -0.704571, -6.93334e-14, -1.86212e-13, 1.23004e-11, -2.4201,
      -1.39111e-13, 1.11167e-12, -6.26898e-11, 4.65755e-12, 1.7833e-13,
      2.6458e-13, 0.141989, 8.09908e-14, -3.50081e-13, 2.30704e-11, -3.92461,
      3.32045e-13, -9.47659e-13, -0.0242444, -0.280382, 0.194573, 2.99591e-11,
      0.309353, 1.11096, -0.810508, -0.666152, -0.102656, -0.931242,
      2.59198e-11, 0.249608, 0.863331, -0.633749, -8.14098e-12, 3.40809e-13,
      -9.9859e-13, 0.252171, 0.11672, 0.199617, 2.93934e-11, 0.807183,
      -0.819287, -0.81491, -0.662062, 0.860059, 0.378789, 2.56054e-11, 0.622636,
      -0.646441, -0.635391, -8.02697e-12, 3.36443e-13, -9.64173e-13, -0.229962,
      0.157553, 0.198382, 2.9661e-11, -1.11338, -0.28529, -0.815998, -0.662891,
      -0.755824, 0.557509, 2.58383e-11, -0.87138, -0.214606, -0.63555,
      -8.08299e-12, -0.000733354, -0.119419, -0.0009883, -0.00296639, 0.287737,
      0.0223445, -0.000165311, -0.000334196, 0.127842, 0.714478, -0.000566861,
      -0.0018141, -0.667062, 0.000398267, 0.00105377, -0.878697, 0.149996,
      -0.016045, 0.00983249, -0.000764951, -0.002296, 0.22271, 0.924462,
      -0.00225927, -0.00456738, 1.74718, -0.15356, 0.000121833, 0.000389896,
      0.302026, -0.000663121, -0.00175455, 1.46304, 0.777358, -0.000733354,
      -0.119419, -0.0225565, -0.269326, -0.0987776, 0.0223445, -0.0322324,
      -0.116014, -0.0429592, -0.236401, 0.0742399, 0.670139, -0.667062,
      0.229906, 0.795514, 0.293958, 0.149996, -0.016045, 0.00983249, -0.0174589,
      -0.20846, -0.0764545, 0.924462, -0.440512, -1.58553, -0.587113, 0.0508086,
      -0.0159561, -0.14403, 0.302026, -0.382798, -1.32454, -0.489444, 0.777358,
      -0.000733354, -0.119419, -0.222339, 0.155973, -0.0950793, 0.0223445,
      0.116767, 0.0302119, -0.0423843, -0.238742, 0.542999, -0.398293,
      -0.667062, -0.804392, -0.199064, 0.292296, 0.149996, -0.016045,
      0.00983249, -0.172092, 0.120724, -0.0735919, 0.924462, 1.59582, 0.412899,
      -0.579255, 0.0513117, -0.116704, 0.0856034, 0.302026, 1.33932, 0.331444,
      -0.486678, 0.777358, -0.000733354, -0.119419, 0.245883, 0.116319,
      -0.0938798, 0.0223445, -0.0843689, 0.0861361, -0.0424982, -0.239336,
      -0.616672, -0.270031, -0.667062, 0.574088, -0.597504, 0.292443, 0.149996,
      -0.016045, 0.00983249, 0.190315, 0.0900316, -0.0726636, 0.924462,
      -1.15305, 1.1772, -0.580812, 0.0514394, 0.132539, 0.0580367, 0.302026,
      -0.955864, 0.994852, -0.486922, 0.777358;
  Eigen::MatrixXd vxc = Eigen::MatrixXd::Zero(17, 17);
  vxc << -1.71426, 0.167004, 3.71196e-14, 5.86163e-14, -1.62646e-12, 0.122215,
      1.13798e-14, -8.69305e-14, 5.83629e-12, 3.18548e-12, -3.11921e-14,
      -4.87405e-14, -0.0904351, -2.79637e-14, 8.08381e-14, -2.52429e-12,
      -0.151944, 0.167004, -0.44154, 4.03427e-14, 4.46379e-14, 5.86056e-13,
      -0.0909547, 1.17684e-14, -3.94129e-15, -3.14496e-12, 7.08364e-13,
      -3.53814e-14, -4.61315e-14, -0.0281286, -1.51545e-14, 6.53366e-14,
      -3.59474e-12, 0.0861955, 3.71231e-14, 4.03601e-14, -0.394774,
      -7.42462e-16, 6.07153e-17, 2.25098e-14, -0.048987, 0.0181965,
      -0.000195258, 0.000513494, 0.123323, -0.0031707, -7.95891e-15,
      0.000614148, -0.000235269, 2.25502e-06, -1.92485e-14, 5.8612e-14,
      4.464e-14, -7.58074e-16, -0.394774, -1.92554e-15, 3.41116e-14, 0.0181975,
      0.0489855, -0.000387129, 0.00158735, 0.00316386, 0.123314, -6.38378e-15,
      -0.00023528, -0.00061412, 5.93669e-06, -4.98282e-14, -1.62646e-12,
      5.8608e-13, 5.0307e-17, -1.92815e-15, -0.394774, -2.75426e-13,
      -4.82304e-05, 0.000430893, 0.052256, -0.123353, 0.000554081, 0.00157365,
      -1.0517e-12, -1.80451e-08, -6.35052e-06, -0.000657643, 1.99361e-12,
      0.122215, -0.0909547, 2.25982e-14, 3.4096e-14, -2.75492e-13, -0.231036,
      -3.85109e-16, -3.85594e-14, 1.51334e-12, -1.29286e-12, -2.1776e-14,
      -3.14152e-14, 0.102952, -2.09277e-14, 5.68677e-14, 9.77943e-12, 0.0984482,
      1.13816e-14, 1.16408e-14, -0.048987, 0.0181975, -4.82304e-05,
      -2.35922e-16, -0.252907, -5.6205e-16, -6.59195e-17, -3.94476e-06,
      0.00243826, -0.000977663, -1.50765e-14, 0.140021, -0.00142595,
      -0.000119316, -1.48631e-14, -8.6965e-14, -3.86827e-15, 0.0181965,
      0.0489855, 0.000430893, -3.86566e-14, -5.6552e-16, -0.252907,
      -5.41581e-15, -1.38337e-05, -0.00097767, -0.00243822, 4.66294e-15,
      0.00142578, 0.140021, -0.000197488, -4.49918e-14, 5.83633e-12,
      -3.14482e-12, -0.000195258, -0.000387129, 0.052256, 1.5134e-12,
      -8.1532e-17, -5.27529e-15, -0.252907, 0.00262693, -1.48708e-06,
      -1.43081e-05, -1.18199e-12, 0.00012132, 0.000196263, 0.140029,
      7.75764e-12, 3.18547e-12, 7.08325e-13, 0.000513494, 0.00158735, -0.123353,
      -1.29288e-12, -3.94476e-06, -1.38337e-05, 0.00262693, -0.311947,
      6.34807e-16, 1.0629e-15, 3.12903e-12, 2.28406e-05, 0.00012763, -0.0331599,
      1.17161e-12, -3.11952e-14, -3.53903e-14, 0.123323, 0.00316386,
      0.000554081, -2.16926e-14, 0.00243826, -0.00097767, -1.48708e-06,
      6.32778e-16, -0.311947, 3.72139e-15, 7.19113e-15, -0.0306508, 0.0126539,
      2.75916e-05, 2.44546e-14, -4.87413e-14, -4.6108e-14, -0.0031707, 0.123314,
      0.00157365, -3.14614e-14, -0.000977663, -0.00243822, -1.43081e-05,
      1.01845e-15, 3.73855e-15, -0.311947, 1.56823e-15, 0.0126539, 0.0306506,
      0.000126688, 3.03234e-14, -0.0904351, -0.0281286, -7.9472e-15,
      -6.45079e-15, -1.05175e-12, 0.102952, -1.50713e-14, 4.60743e-15,
      -1.18208e-12, 3.12903e-12, 7.26415e-15, 1.63389e-15, -0.406695,
      1.4086e-15, -2.68275e-14, 6.20841e-12, -0.03473, -2.7965e-14,
      -1.52253e-14, 0.000614148, -0.00023528, -1.80451e-08, -2.08271e-14,
      0.140021, 0.00142578, 0.00012132, 2.28406e-05, -0.0306508, 0.0126539,
      1.4615e-15, -0.410651, -6.10623e-16, 4.51028e-17, 1.39472e-15,
      8.08144e-14, 6.54684e-14, -0.000235269, -0.00061412, -6.35052e-06,
      5.67532e-14, -0.00142595, 0.140021, 0.000196263, 0.00012763, 0.0126539,
      0.0306506, -2.68154e-14, -6.07153e-16, -0.410651, -8.51749e-15,
      -3.30638e-14, -2.52428e-12, -3.59461e-12, 2.25502e-06, 5.93669e-06,
      -0.000657643, 9.77951e-12, -0.000119316, -0.000197488, 0.140029,
      -0.0331599, 2.75916e-05, 0.000126688, 6.20841e-12, 5.55112e-17,
      -8.43423e-15, -0.410651, 4.25785e-12, -0.151944, 0.0861955, -1.92084e-14,
      -4.98197e-14, 1.99354e-12, 0.0984482, -1.46012e-14, -4.48374e-14,
      7.75783e-12, 1.17156e-12, 2.4393e-14, 3.03816e-14, -0.03473, 1.55258e-15,
      -3.28227e-14, 4.258e-12, -0.352989;

  Orbitals orbitals;
  orbitals.QMAtoms().LoadFromFile("molecule.xyz");
  BasisSet basis;
  basis.Load("3-21G.xml");
  orbitals.setDFTbasisName("3-21G.xml");
  AOBasis aobasis;
  aobasis.Fill(basis, orbitals.QMAtoms());
  orbitals.setBasisSetSize(17);
  orbitals.setNumberOfOccupiedLevels(4);
  orbitals.MOs().eigenvectors() = mo_eigenvectors;
  Logger log;
  TCMatrix_gwbse Mmn{log};
  Mmn.Initialize(aobasis.AOBasisSize(), 0, 16, 0, 16);
  Mmn.Fill(aobasis, aobasis, mo_eigenvectors);

  GW::options opt;
  opt.ScaHFX = 0;
  opt.homo = 4;
  opt.qpmax = 16;
  opt.qpmin = 0;
  opt.rpamax = 16;
  opt.rpamin = 0;
  opt.gw_sc_max_iterations = 1;
  opt.eta = 1e-3;
  GW gw(log, Mmn, vxc, mo_eigenvalues);
  gw.configure(opt);

  Eigen::MatrixXd ref = Eigen::MatrixXd::Zero(17, 17);
  ref << -11.4075, -0.000332076, -2.98562e-07, -1.444e-07, 3.56698e-08,
      0.00415282, -2.59416e-07, 3.42809e-07, 2.37553e-08, -1.12992e-07,
      1.90622e-07, 7.34527e-08, -0.00225178, -9.78891e-07, 2.05282e-07,
      -8.26113e-07, -0.0537384, -0.000332076, -1.01045, -9.0233e-07,
      -5.36336e-07, 2.09522e-07, 0.00201438, -4.15248e-07, 6.11818e-07,
      1.77817e-07, -4.04586e-07, 4.29371e-07, 3.58778e-07, 0.00539939,
      -1.90892e-06, 4.37732e-07, -1.69582e-06, 0.00539329, -2.98562e-07,
      -9.0233e-07, -0.620146, 1.97176e-07, -7.69219e-07, 1.12519e-07, 0.0208017,
      -0.00772665, 8.27053e-05, 8.52774e-05, 0.0206388, -0.000530697,
      -3.25714e-07, 0.00634783, -0.00243228, 2.30183e-05, -1.20011e-07,
      -1.444e-07, -5.36336e-07, 1.97176e-07, -0.620147, -6.94777e-07,
      3.49514e-08, -0.00772665, -0.0208006, 0.000164225, 0.000265677,
      0.000529342, 0.020638, -1.52421e-07, -0.00243197, -0.00634853,
      6.17472e-05, 9.44e-08, 3.56698e-08, 2.09522e-07, -7.69219e-07,
      -6.94777e-07, -0.620149, -5.8171e-08, 2.06511e-05, -0.000183069,
      -0.0221897, -0.0206448, 9.304e-05, 0.000263363, 9.48884e-08, -3.68494e-07,
      -6.53318e-05, -0.00679866, 1.96631e-07, 0.00415282, 0.00201438,
      1.12519e-07, 3.49514e-08, -5.8171e-08, 0.261183, -7.74632e-08, 4.7791e-08,
      -3.39764e-08, -1.29803e-08, -1.37184e-08, -8.65211e-09, 0.0282484,
      -1.08964e-06, 2.33312e-07, -1.06382e-06, -0.0110563, -2.59416e-07,
      -4.15248e-07, 0.0208017, -0.00772665, 2.06511e-05, -7.74632e-08, 0.348343,
      -4.66202e-08, -1.8013e-08, -1.68356e-05, 0.0104315, -0.00418275,
      9.38435e-08, -0.0337371, 0.000343173, 2.90984e-05, -2.15296e-09,
      3.42809e-07, 6.11818e-07, -0.00772665, -0.0208006, -0.000183069,
      4.7791e-08, -4.66202e-08, 0.348343, 6.05195e-08, -5.92169e-05,
      -0.00418276, -0.0104315, 7.04968e-08, -0.000343952, -0.0337369,
      4.7414e-05, -3.6549e-08, 2.37553e-08, 1.77817e-07, 8.27053e-05,
      0.000164225, -0.0221897, -3.39764e-08, -1.8013e-08, 6.05195e-08, 0.348342,
      0.0112388, -6.25695e-06, -6.11456e-05, 8.74033e-08, -2.86503e-05,
      -4.75377e-05, -0.0337397, -4.59301e-08, -1.12992e-07, -4.04586e-07,
      8.52774e-05, 0.000265677, -0.0206448, -1.29803e-08, -1.68356e-05,
      -5.92169e-05, 0.0112388, 0.920796, 1.25585e-07, 3.1022e-08, -4.1467e-08,
      -1.39666e-05, -7.56229e-05, 0.0197026, 1.82784e-07, 1.90622e-07,
      4.29371e-07, 0.0206388, 0.000529342, 9.304e-05, -1.37184e-08, 0.0104315,
      -0.00418276, -6.25695e-06, 1.25585e-07, 0.920796, 5.39747e-08,
      3.64479e-07, 0.0182118, -0.00751815, -1.66449e-05, -6.84045e-08,
      7.34527e-08, 3.58778e-07, -0.000530697, 0.020638, 0.000263363,
      -8.65211e-09, -0.00418275, -0.0104315, -6.11456e-05, 3.1022e-08,
      5.39747e-08, 0.920795, 2.06775e-07, -0.00751882, -0.0182113, -7.54346e-05,
      -1.30207e-07, -0.00225178, 0.00539939, -3.25714e-07, -1.52421e-07,
      9.48884e-08, 0.0282484, 9.38435e-08, 7.04968e-08, 8.74033e-08,
      -4.1467e-08, 3.64479e-07, 2.06775e-07, 1.13224, 3.48068e-07, -8.5456e-08,
      3.72352e-07, -0.0624424, -9.78891e-07, -1.90892e-06, 0.00634783,
      -0.00243197, -3.68494e-07, -1.08964e-06, -0.0337371, -0.000343952,
      -2.86503e-05, -1.39666e-05, 0.0182118, -0.00751882, 3.48068e-07, 1.29628,
      1.7215e-07, -1.96562e-07, -3.24891e-07, 2.05282e-07, 4.37732e-07,
      -0.00243228, -0.00634853, -6.53318e-05, 2.33312e-07, 0.000343173,
      -0.0337369, -4.75377e-05, -7.56229e-05, -0.00751815, -0.0182113,
      -8.5456e-08, 1.7215e-07, 1.29628, 1.85216e-08, 2.47821e-07, -8.26113e-07,
      -1.69582e-06, 2.30183e-05, 6.17472e-05, -0.00679866, -1.06382e-06,
      2.90984e-05, 4.7414e-05, -0.0337397, 0.0197026, -1.66449e-05,
      -7.54346e-05, 3.72352e-07, -1.96562e-07, 1.85216e-08, 1.29628,
      -1.65205e-07, -0.0537384, 0.00539329, -1.20011e-07, 9.44e-08, 1.96631e-07,
      -0.0110563, -2.15296e-09, -3.6549e-08, -4.59301e-08, 1.82784e-07,
      -6.84045e-08, -1.30207e-07, -0.0624424, -3.24891e-07, 2.47821e-07,
      -1.65205e-07, 1.96234;

  gw.CalculateGWPerturbation();
  Eigen::VectorXd diag = gw.getGWAResults();
  bool check_diag = ref.diagonal().isApprox(diag, 1e-4);
  if (!check_diag) {
    cout << "GW energies" << endl;
    cout << diag << endl;
    cout << "GW energies ref" << endl;
    cout << ref.diagonal() << endl;
  }
  BOOST_CHECK_EQUAL(check_diag, true);

  gw.CalculateHQP();
  Eigen::MatrixXd offdiag = gw.getHQP();

  bool check_offdiag = ref.isApprox(offdiag, 1e-4);
  if (!check_offdiag) {
    cout << "GW energies" << endl;
    cout << offdiag << endl;
    cout << "GW energies ref" << endl;
    cout << ref << endl;
  }
  BOOST_CHECK_EQUAL(check_offdiag, true);
}

BOOST_AUTO_TEST_CASE(gw_full_QP_grid) {

  ofstream xyzfile("molecule.xyz");
  xyzfile << " 5" << endl;
  xyzfile << " methane" << endl;
  xyzfile << " C            .000000     .000000     .000000" << endl;
  xyzfile << " H            .629118     .629118     .629118" << endl;
  xyzfile << " H           -.629118    -.629118     .629118" << endl;
  xyzfile << " H            .629118    -.629118    -.629118" << endl;
  xyzfile << " H           -.629118     .629118    -.629118" << endl;
  xyzfile.close();

  ofstream basisfile("3-21G.xml");
  basisfile << "<basis name=\"3-21G\">" << endl;
  basisfile << "  <element name=\"H\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"5.447178e+00\">" << endl;
  basisfile << "        <contractions factor=\"1.562850e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"8.245470e-01\">" << endl;
  basisfile << "        <contractions factor=\"9.046910e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"1.831920e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "  <element name=\"C\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"1.722560e+02\">" << endl;
  basisfile << "        <contractions factor=\"6.176690e-02\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"2.591090e+01\">" << endl;
  basisfile << "        <contractions factor=\"3.587940e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"5.533350e+00\">" << endl;
  basisfile << "        <contractions factor=\"7.007130e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << endl;
  basisfile << "      <constant decay=\"3.664980e+00\">" << endl;
  basisfile << "        <contractions factor=\"-3.958970e-01\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"2.364600e-01\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"7.705450e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.215840e+00\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"8.606190e-01\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << endl;
  basisfile << "      <constant decay=\"1.958570e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "</basis>" << endl;
  basisfile.close();

  Eigen::VectorXd mo_eigenvalues = Eigen::VectorXd::Zero(17);
  mo_eigenvalues << -10.6784, -0.746424, -0.394948, -0.394948, -0.394948,
      0.165212, 0.227713, 0.227713, 0.227713, 0.763971, 0.763971, 0.763971,
      1.05054, 1.13372, 1.13372, 1.13372, 1.72964;
  Eigen::MatrixXd mo_eigenvectors = Eigen::MatrixXd::Zero(17, 17);
  mo_eigenvectors << -0.990744, 0.205494, 1.12419e-14, 1.45929e-14, 3.75871e-15,
      0.169799, 8.81326e-15, -7.59696e-14, 4.96934e-12, 3.56824e-13,
      -4.00071e-15, -9.34192e-15, -0.116538, -1.31882e-14, 3.3298e-14,
      -4.67337e-12, 0.0231349, -0.089086, -0.203519, 1.36072e-14, 4.0936e-14,
      -2.48734e-12, -0.0756531, 8.87485e-15, 1.17892e-14, -6.66069e-12,
      -3.61695e-12, -5.11501e-14, -5.65398e-14, 0.160728, 4.94743e-15,
      8.66668e-14, -1.06301e-11, 1.99744, -1.87806e-15, -2.96718e-13,
      -0.0253717, -0.29342, 0.203621, 1.37293e-12, 0.0616491, 0.221397,
      -0.161522, 0.405432, 0.0624782, 0.566771, -2.44075e-11, -0.205261,
      -0.709949, 0.521155, 5.26171e-12, -1.61071e-14, -4.43669e-13, 0.263897,
      0.122148, 0.208899, 1.16501e-12, 0.160859, -0.163271, -0.162399, 0.402943,
      -0.523447, -0.230538, -2.4163e-11, -0.512016, 0.531592, 0.522505,
      5.1261e-12, -8.74496e-15, -3.64573e-13, -0.240655, 0.164879, 0.207606,
      1.26241e-12, -0.221879, -0.0568539, -0.162616, 0.403447, 0.460008,
      -0.33931, -2.43511e-11, 0.716567, 0.176478, 0.522636, 5.19033e-12,
      0.0813118, -0.704571, -6.93334e-14, -1.86212e-13, 1.23004e-11, -2.4201,
      -1.39111e-13, 1.11167e-12, -6.26898e-11, 4.65755e-12, 1.7833e-13,
      2.6458e-13, 0.141989, 8.09908e-14, -3.50081e-13, 2.30704e-11, -3.92461,
      3.32045e-13, -9.47659e-13, -0.0242444, -0.280382, 0.194573, 2.99591e-11,
      0.309353, 1.11096, -0.810508, -0.666152, -0.102656, -0.931242,
      2.59198e-11, 0.249608, 0.863331, -0.633749, -8.14098e-12, 3.40809e-13,
      -9.9859e-13, 0.252171, 0.11672, 0.199617, 2.93934e-11, 0.807183,
      -0.819287, -0.81491, -0.662062, 0.860059, 0.378789, 2.56054e-11, 0.622636,
      -0.646441, -0.635391, -8.02697e-12, 3.36443e-13, -9.64173e-13, -0.229962,
      0.157553, 0.198382, 2.9661e-11, -1.11338, -0.28529, -0.815998, -0.662891,
      -0.755824, 0.557509, 2.58383e-11, -0.87138, -0.214606, -0.63555,
      -8.08299e-12, -0.000733354, -0.119419, -0.0009883, -0.00296639, 0.287737,
      0.0223445, -0.000165311, -0.000334196, 0.127842, 0.714478, -0.000566861,
      -0.0018141, -0.667062, 0.000398267, 0.00105377, -0.878697, 0.149996,
      -0.016045, 0.00983249, -0.000764951, -0.002296, 0.22271, 0.924462,
      -0.00225927, -0.00456738, 1.74718, -0.15356, 0.000121833, 0.000389896,
      0.302026, -0.000663121, -0.00175455, 1.46304, 0.777358, -0.000733354,
      -0.119419, -0.0225565, -0.269326, -0.0987776, 0.0223445, -0.0322324,
      -0.116014, -0.0429592, -0.236401, 0.0742399, 0.670139, -0.667062,
      0.229906, 0.795514, 0.293958, 0.149996, -0.016045, 0.00983249, -0.0174589,
      -0.20846, -0.0764545, 0.924462, -0.440512, -1.58553, -0.587113, 0.0508086,
      -0.0159561, -0.14403, 0.302026, -0.382798, -1.32454, -0.489444, 0.777358,
      -0.000733354, -0.119419, -0.222339, 0.155973, -0.0950793, 0.0223445,
      0.116767, 0.0302119, -0.0423843, -0.238742, 0.542999, -0.398293,
      -0.667062, -0.804392, -0.199064, 0.292296, 0.149996, -0.016045,
      0.00983249, -0.172092, 0.120724, -0.0735919, 0.924462, 1.59582, 0.412899,
      -0.579255, 0.0513117, -0.116704, 0.0856034, 0.302026, 1.33932, 0.331444,
      -0.486678, 0.777358, -0.000733354, -0.119419, 0.245883, 0.116319,
      -0.0938798, 0.0223445, -0.0843689, 0.0861361, -0.0424982, -0.239336,
      -0.616672, -0.270031, -0.667062, 0.574088, -0.597504, 0.292443, 0.149996,
      -0.016045, 0.00983249, 0.190315, 0.0900316, -0.0726636, 0.924462,
      -1.15305, 1.1772, -0.580812, 0.0514394, 0.132539, 0.0580367, 0.302026,
      -0.955864, 0.994852, -0.486922, 0.777358;
  Eigen::MatrixXd vxc = Eigen::MatrixXd::Zero(17, 17);
  vxc << -1.71426, 0.167004, 3.71196e-14, 5.86163e-14, -1.62646e-12, 0.122215,
      1.13798e-14, -8.69305e-14, 5.83629e-12, 3.18548e-12, -3.11921e-14,
      -4.87405e-14, -0.0904351, -2.79637e-14, 8.08381e-14, -2.52429e-12,
      -0.151944, 0.167004, -0.44154, 4.03427e-14, 4.46379e-14, 5.86056e-13,
      -0.0909547, 1.17684e-14, -3.94129e-15, -3.14496e-12, 7.08364e-13,
      -3.53814e-14, -4.61315e-14, -0.0281286, -1.51545e-14, 6.53366e-14,
      -3.59474e-12, 0.0861955, 3.71231e-14, 4.03601e-14, -0.394774,
      -7.42462e-16, 6.07153e-17, 2.25098e-14, -0.048987, 0.0181965,
      -0.000195258, 0.000513494, 0.123323, -0.0031707, -7.95891e-15,
      0.000614148, -0.000235269, 2.25502e-06, -1.92485e-14, 5.8612e-14,
      4.464e-14, -7.58074e-16, -0.394774, -1.92554e-15, 3.41116e-14, 0.0181975,
      0.0489855, -0.000387129, 0.00158735, 0.00316386, 0.123314, -6.38378e-15,
      -0.00023528, -0.00061412, 5.93669e-06, -4.98282e-14, -1.62646e-12,
      5.8608e-13, 5.0307e-17, -1.92815e-15, -0.394774, -2.75426e-13,
      -4.82304e-05, 0.000430893, 0.052256, -0.123353, 0.000554081, 0.00157365,
      -1.0517e-12, -1.80451e-08, -6.35052e-06, -0.000657643, 1.99361e-12,
      0.122215, -0.0909547, 2.25982e-14, 3.4096e-14, -2.75492e-13, -0.231036,
      -3.85109e-16, -3.85594e-14, 1.51334e-12, -1.29286e-12, -2.1776e-14,
      -3.14152e-14, 0.102952, -2.09277e-14, 5.68677e-14, 9.77943e-12, 0.0984482,
      1.13816e-14, 1.16408e-14, -0.048987, 0.0181975, -4.82304e-05,
      -2.35922e-16, -0.252907, -5.6205e-16, -6.59195e-17, -3.94476e-06,
      0.00243826, -0.000977663, -1.50765e-14, 0.140021, -0.00142595,
      -0.000119316, -1.48631e-14, -8.6965e-14, -3.86827e-15, 0.0181965,
      0.0489855, 0.000430893, -3.86566e-14, -5.6552e-16, -0.252907,
      -5.41581e-15, -1.38337e-05, -0.00097767, -0.00243822, 4.66294e-15,
      0.00142578, 0.140021, -0.000197488, -4.49918e-14, 5.83633e-12,
      -3.14482e-12, -0.000195258, -0.000387129, 0.052256, 1.5134e-12,
      -8.1532e-17, -5.27529e-15, -0.252907, 0.00262693, -1.48708e-06,
      -1.43081e-05, -1.18199e-12, 0.00012132, 0.000196263, 0.140029,
      7.75764e-12, 3.18547e-12, 7.08325e-13, 0.000513494, 0.00158735, -0.123353,
      -1.29288e-12, -3.94476e-06, -1.38337e-05, 0.00262693, -0.311947,
      6.34807e-16, 1.0629e-15, 3.12903e-12, 2.28406e-05, 0.00012763, -0.0331599,
      1.17161e-12, -3.11952e-14, -3.53903e-14, 0.123323, 0.00316386,
      0.000554081, -2.16926e-14, 0.00243826, -0.00097767, -1.48708e-06,
      6.32778e-16, -0.311947, 3.72139e-15, 7.19113e-15, -0.0306508, 0.0126539,
      2.75916e-05, 2.44546e-14, -4.87413e-14, -4.6108e-14, -0.0031707, 0.123314,
      0.00157365, -3.14614e-14, -0.000977663, -0.00243822, -1.43081e-05,
      1.01845e-15, 3.73855e-15, -0.311947, 1.56823e-15, 0.0126539, 0.0306506,
      0.000126688, 3.03234e-14, -0.0904351, -0.0281286, -7.9472e-15,
      -6.45079e-15, -1.05175e-12, 0.102952, -1.50713e-14, 4.60743e-15,
      -1.18208e-12, 3.12903e-12, 7.26415e-15, 1.63389e-15, -0.406695,
      1.4086e-15, -2.68275e-14, 6.20841e-12, -0.03473, -2.7965e-14,
      -1.52253e-14, 0.000614148, -0.00023528, -1.80451e-08, -2.08271e-14,
      0.140021, 0.00142578, 0.00012132, 2.28406e-05, -0.0306508, 0.0126539,
      1.4615e-15, -0.410651, -6.10623e-16, 4.51028e-17, 1.39472e-15,
      8.08144e-14, 6.54684e-14, -0.000235269, -0.00061412, -6.35052e-06,
      5.67532e-14, -0.00142595, 0.140021, 0.000196263, 0.00012763, 0.0126539,
      0.0306506, -2.68154e-14, -6.07153e-16, -0.410651, -8.51749e-15,
      -3.30638e-14, -2.52428e-12, -3.59461e-12, 2.25502e-06, 5.93669e-06,
      -0.000657643, 9.77951e-12, -0.000119316, -0.000197488, 0.140029,
      -0.0331599, 2.75916e-05, 0.000126688, 6.20841e-12, 5.55112e-17,
      -8.43423e-15, -0.410651, 4.25785e-12, -0.151944, 0.0861955, -1.92084e-14,
      -4.98197e-14, 1.99354e-12, 0.0984482, -1.46012e-14, -4.48374e-14,
      7.75783e-12, 1.17156e-12, 2.4393e-14, 3.03816e-14, -0.03473, 1.55258e-15,
      -3.28227e-14, 4.258e-12, -0.352989;

  Orbitals orbitals;
  orbitals.QMAtoms().LoadFromFile("molecule.xyz");
  BasisSet basis;
  basis.Load("3-21G.xml");
  orbitals.setDFTbasisName("3-21G.xml");
  AOBasis aobasis;
  aobasis.Fill(basis, orbitals.QMAtoms());
  orbitals.setBasisSetSize(17);
  orbitals.setNumberOfOccupiedLevels(4);
  orbitals.MOs().eigenvectors() = mo_eigenvectors;
  Logger log;
  TCMatrix_gwbse Mmn{log};
  Mmn.Initialize(aobasis.AOBasisSize(), 0, 16, 0, 16);
  Mmn.Fill(aobasis, aobasis, mo_eigenvectors);

  GW::options opt;
  opt.ScaHFX = 0;
  opt.homo = 4;
  opt.qpmax = 16;
  opt.qpmin = 0;
  opt.rpamax = 16;
  opt.rpamin = 0;
  opt.gw_sc_max_iterations = 1;
  opt.qp_solver = "grid";
  opt.eta = 1e-3;
  GW gw(log, Mmn, vxc, mo_eigenvalues);
  gw.configure(opt);

  gw.CalculateGWPerturbation();
  Eigen::VectorXd diag = gw.getGWAResults();

  Eigen::MatrixXd ref = Eigen::MatrixXd::Zero(17, 17);
  ref << -11.4075, -0.000331803, -2.98562e-07, -1.44399e-07, 3.56698e-08,
      0.004153, -2.59416e-07, 3.42809e-07, 2.37555e-08, -1.12993e-07,
      1.90621e-07, 7.34521e-08, -0.00225169, -9.78896e-07, 2.05284e-07,
      -8.26116e-07, -0.0537394, -0.000331803, -1.01045, -9.0233e-07,
      -5.36336e-07, 2.09522e-07, 0.00201439, -4.15249e-07, 6.11818e-07,
      1.77817e-07, -4.04586e-07, 4.29371e-07, 3.58778e-07, 0.00539987,
      -1.9089e-06, 4.37728e-07, -1.6958e-06, 0.00539247, -2.98562e-07,
      -9.0233e-07, -0.620145, 1.97176e-07, -7.69219e-07, 1.12519e-07, 0.0208017,
      -0.00772665, 8.27053e-05, 8.52774e-05, 0.0206388, -0.000530697,
      -3.25724e-07, 0.00634825, -0.00243244, 2.30204e-05, -1.20016e-07,
      -1.44399e-07, -5.36336e-07, 1.97176e-07, -0.620147, -6.94777e-07,
      3.49514e-08, -0.00772665, -0.0208006, 0.000164225, 0.000265677,
      0.000529343, 0.020638, -1.5243e-07, -0.00243213, -0.00634896, 6.17529e-05,
      9.43974e-08, 3.56698e-08, 2.09522e-07, -7.69219e-07, -6.94777e-07,
      -0.620148, -5.8171e-08, 2.06511e-05, -0.000183069, -0.0221897, -0.0206448,
      9.30401e-05, 0.000263363, 9.48807e-08, -3.68484e-07, -6.53363e-05,
      -0.0067993, 1.96629e-07, 0.004153, 0.00201439, 1.12519e-07, 3.49514e-08,
      -5.8171e-08, 0.261183, -7.74631e-08, 4.7791e-08, -3.39764e-08,
      -1.29803e-08, -1.37184e-08, -8.65211e-09, 0.0282494, -1.08969e-06,
      2.33331e-07, -1.06388e-06, -0.0110577, -2.59416e-07, -4.15249e-07,
      0.0208017, -0.00772665, 2.06511e-05, -7.74631e-08, 0.348343, -4.66202e-08,
      -1.8013e-08, -1.68356e-05, 0.0104315, -0.00418275, 9.38381e-08,
      -0.0337379, 0.000343182, 2.90994e-05, -2.1529e-09, 3.42809e-07,
      6.11818e-07, -0.00772665, -0.0208006, -0.000183069, 4.7791e-08,
      -4.66202e-08, 0.348343, 6.05195e-08, -5.92169e-05, -0.00418276,
      -0.0104315, 7.05055e-08, -0.000343961, -0.0337378, 4.74156e-05,
      -3.65502e-08, 2.37555e-08, 1.77817e-07, 8.27053e-05, 0.000164225,
      -0.0221897, -3.39764e-08, -1.8013e-08, 6.05195e-08, 0.348343, 0.0112388,
      -6.25696e-06, -6.11457e-05, 8.74054e-08, -2.8651e-05, -4.75389e-05,
      -0.0337408, -4.59305e-08, -1.12993e-07, -4.04586e-07, 8.52774e-05,
      0.000265677, -0.0206448, -1.29803e-08, -1.68356e-05, -5.92169e-05,
      0.0112388, 0.920796, 1.25585e-07, 3.1022e-08, -4.14705e-08, -1.39666e-05,
      -7.56233e-05, 0.0197027, 1.82784e-07, 1.90621e-07, 4.29371e-07, 0.0206388,
      0.000529343, 9.30401e-05, -1.37184e-08, 0.0104315, -0.00418276,
      -6.25696e-06, 1.25585e-07, 0.920796, 5.39746e-08, 3.64492e-07, 0.0182119,
      -0.00751818, -1.6645e-05, -6.83932e-08, 7.34521e-08, 3.58778e-07,
      -0.000530697, 0.020638, 0.000263363, -8.65211e-09, -0.00418275,
      -0.0104315, -6.11457e-05, 3.1022e-08, 5.39746e-08, 0.920796, 2.06784e-07,
      -0.00751885, -0.0182114, -7.54351e-05, -1.30202e-07, -0.00225169,
      0.00539987, -3.25724e-07, -1.5243e-07, 9.48807e-08, 0.0282494,
      9.38381e-08, 7.05055e-08, 8.74054e-08, -4.14705e-08, 3.64492e-07,
      2.06784e-07, 1.13224, 3.48114e-07, -8.54848e-08, 3.72384e-07, -0.0624423,
      -9.78896e-07, -1.9089e-06, 0.00634825, -0.00243213, -3.68484e-07,
      -1.08969e-06, -0.0337379, -0.000343961, -2.8651e-05, -1.39666e-05,
      0.0182119, -0.00751885, 3.48114e-07, 1.29628, 1.72158e-07, -1.96578e-07,
      -3.24921e-07, 2.05284e-07, 4.37728e-07, -0.00243244, -0.00634896,
      -6.53363e-05, 2.33331e-07, 0.000343182, -0.0337378, -4.75389e-05,
      -7.56233e-05, -0.00751818, -0.0182114, -8.54848e-08, 1.72158e-07, 1.29628,
      1.85231e-08, 2.47837e-07, -8.26116e-07, -1.6958e-06, 2.30204e-05,
      6.17529e-05, -0.0067993, -1.06388e-06, 2.90994e-05, 4.74156e-05,
      -0.0337408, 0.0197027, -1.6645e-05, -7.54351e-05, 3.72384e-07,
      -1.96578e-07, 1.85231e-08, 1.29628, -1.65237e-07, -0.0537394, 0.00539247,
      -1.20016e-07, 9.43974e-08, 1.96629e-07, -0.0110577, -2.1529e-09,
      -3.65502e-08, -4.59305e-08, 1.82784e-07, -6.83932e-08, -1.30202e-07,
      -0.0624423, -3.24921e-07, 2.47837e-07, -1.65237e-07, 1.96233;

  bool check_diag = ref.diagonal().isApprox(diag, 1e-4);
  if (!check_diag) {
    cout << "GW energies" << endl;
    cout << diag << endl;
    cout << "GW energies ref" << endl;
    cout << ref.diagonal() << endl;
  }
  BOOST_CHECK_EQUAL(check_diag, true);

  gw.CalculateHQP();
  Eigen::MatrixXd offdiag = gw.getHQP();

  bool check_offdiag = ref.isApprox(offdiag, 1e-4);
  if (!check_offdiag) {
    cout << "GW energies" << endl;
    cout << offdiag << endl;
    cout << "GW energies ref" << endl;
    cout << ref << endl;
  }
  BOOST_CHECK_EQUAL(check_offdiag, true);
}

BOOST_AUTO_TEST_SUITE_END()
